##################################################################
## Replicating Burmila, "The Electoral College"                 ##
## Krista Loose and Vanessa Williamson                          ##
##################################################################

##################################################################
## Install and load packages needed for following code to run.  ##
##################################################################

##install.packages("plyr")			
##install.packages("reshape")
##install.pacakges("foreign")
##install.packages("maps")
##install.packages("xtable")

library(plyr)
library(reshape)
library(foreign)
library(maps)
library(xtable)

##################################################################
## Download data.                                               ##
##################################################################

A <- read.dta("votes.dta")
B <- read.dta("census.dta")

##################################################################
## Calculate actual democratic share of 2-party vote in each    ##
## election.                                                    ##
## 												                ##
## Data from 1996-2008 are available on the FEC website at      ##
## http://www.fec.gov/pubrec/electionresults.shtml              ##  
## Data from 1948-1992 were obtained from		                ##
## www.uselectionatlas.com and will be verified with government ##
## records to ensure accuracy.  				                ##
## Note that data are not available for all states for all      ##
## years; Alaska and Hawaii did not become states until 1959    ##
## and DC did not obtain electoral votes until the 1964         ##
## election.                                                    ##  
##################################################################

attach(A)
A$ave2008 <- d2008/(d2008+r2008)
A$ave2004 <- d2004/(d2004+r2004) 
A$ave2000 <- d2000/(d2000+r2000) 
A$ave1996 <- d1996/(d1996+r1996)
A$ave1992 <- d1992/(d1992+r1992)
A$ave1988 <- d1988/(d1988+r1988)
A$ave1984 <- d1984/(d1984+r1984)
A$ave1980 <- d1980/(d1980+r1980)
A$ave1976 <- d1976/(d1976+r1976)
A$ave1972 <- d1972/(d1972+r1972)
A$ave1968 <- d1968/(d1968+r1968)
A$ave1964 <- d1964/(d1964+r1964)
A$ave1960 <- d1960/(d1960+r1960)
A$ave1956 <- d1956/(d1956+r1956)
A$ave1952 <- d1952/(d1952+r1952)
A$ave1948 <- d1948/(d1948+r1948)
A$ave1944 <- d1944/(d1944+r1944)
A$ave1940 <- d1940/(d1940+r1940)
A$ave1936 <- d1936/(d1936+r1936)
A$ave1932 <- d1932/(d1932+r1932)
detach(A)

##################################################################
## Calculate normal vote predictions for elections from 1968 to ##
## 2008 using Burmila's WNV equation for previous 5 elections   ## 
## (equivalent to Burmila's predictions for 2012).  These are   ##
## called p.1968 - p.2008. Calculate normal vote predictions    ##
## for elections from 1972 to 2008 using Burmila's WNV equation ##
## for elections starting 6 elections previously and ending 2   ##
## elections previously (equivalent to Burmila's predictions    ##
## for 2016).  These are called p4.1972 - p4.2008.  Continue    ##
## this pattern until you have predictions equivalent to        ##
## Burmila's prediction for 2028.  These are called p16.1984 -  ##
## p16.2008 and use election data from 20 years earlier.        ##
##################################################################

attach(A)

A$p.2008 <- (ave2004*5 + ave2000*4 + ave1996*3 + ave1992*2 + ave1988)/15  ## equivalent to 2012
A$p.2004 <- (ave2000*5 + ave1996*4 + ave1992*3 + ave1988*2 + ave1984)/15  
A$p.2000 <- (ave1996*5 + ave1992*4 + ave1988*3 + ave1984*2 + ave1980)/15 
A$p.1996 <- (ave1992*5 + ave1988*4 + ave1984*3 + ave1980*2 + ave1976)/15
A$p.1992 <- (ave1988*5 + ave1984*4 + ave1980*3 + ave1976*2 + ave1972)/15
A$p.1988 <- (ave1984*5 + ave1980*4 + ave1976*3 + ave1972*2 + ave1968)/15
A$p.1984 <- (ave1980*5 + ave1976*4 + ave1972*3 + ave1968*2 + ave1964)/15
A$p.1980 <- (ave1976*5 + ave1972*4 + ave1968*3 + ave1964*2 + ave1960)/15
A$p.1976 <- (ave1972*5 + ave1968*4 + ave1964*3 + ave1960*2 + ave1956)/15
A$p.1972 <- (ave1968*5 + ave1964*4 + ave1960*3 + ave1956*2 + ave1952)/15
A$p.1968 <- (ave1964*5 + ave1960*4 + ave1956*3 + ave1952*2 + ave1948)/15
A$p.1964 <- (ave1960*5 + ave1956*4 + ave1952*3 + ave1948*2 + ave1944)/15
A$p.1960 <- (ave1956*5 + ave1952*4 + ave1948*3 + ave1944*2 + ave1940)/15
A$p.1956 <- (ave1952*5 + ave1948*4 + ave1944*3 + ave1940*2 + ave1936)/15
A$p.1952 <- (ave1948*5 + ave1944*4 + ave1940*3 + ave1936*2 + ave1932)/15

A$p4.2008 <- (ave2000*5 + ave1996*4 + ave1992*3 + ave1988*2 + ave1984)/15  ## equivalent to 2016
A$p4.2004 <- (ave1996*5 + ave1992*4 + ave1988*3 + ave1984*2 + ave1980)/15  
A$p4.2000 <- (ave1992*5 + ave1988*4 + ave1984*3 + ave1980*2 + ave1976)/15
A$p4.1996 <- (ave1988*5 + ave1984*4 + ave1980*3 + ave1976*2 + ave1972)/15
A$p4.1992 <- (ave1984*5 + ave1980*4 + ave1976*3 + ave1972*2 + ave1968)/15
A$p4.1988 <- (ave1980*5 + ave1976*4 + ave1972*3 + ave1968*2 + ave1964)/15
A$p4.1984 <- (ave1976*5 + ave1972*4 + ave1968*3 + ave1964*2 + ave1960)/15
A$p4.1980 <- (ave1972*5 + ave1968*4 + ave1964*3 + ave1960*2 + ave1956)/15
A$p4.1976 <- (ave1968*5 + ave1964*4 + ave1960*3 + ave1956*2 + ave1952)/15
A$p4.1972 <- (ave1964*5 + ave1960*4 + ave1956*3 + ave1952*2 + ave1948)/15
A$p4.1968 <- (ave1960*5 + ave1956*4 + ave1952*3 + ave1948*2 + ave1944)/15
A$p4.1964 <- (ave1956*5 + ave1952*4 + ave1948*3 + ave1944*2 + ave1940)/15
A$p4.1960 <- (ave1952*5 + ave1948*4 + ave1944*3 + ave1940*2 + ave1936)/15
A$p4.1956 <- (ave1948*5 + ave1944*4 + ave1940*3 + ave1936*2 + ave1932)/15

A$p8.2008 <- (ave1996*5 + ave1992*4 + ave1988*3 + ave1984*2 + ave1980)/15  ## equivalent to 2020
A$p8.2004 <- (ave1992*5 + ave1988*4 + ave1984*3 + ave1980*2 + ave1976)/15
A$p8.2000 <- (ave1988*5 + ave1984*4 + ave1980*3 + ave1976*2 + ave1972)/15
A$p8.1996 <- (ave1984*5 + ave1980*4 + ave1976*3 + ave1972*2 + ave1968)/15
A$p8.1992 <- (ave1980*5 + ave1976*4 + ave1972*3 + ave1968*2 + ave1964)/15
A$p8.1988 <- (ave1976*5 + ave1972*4 + ave1968*3 + ave1964*2 + ave1960)/15
A$p8.1984 <- (ave1972*5 + ave1968*4 + ave1964*3 + ave1960*2 + ave1956)/15
A$p8.1980 <- (ave1968*5 + ave1964*4 + ave1960*3 + ave1956*2 + ave1952)/15
A$p8.1976 <- (ave1964*5 + ave1960*4 + ave1956*3 + ave1952*2 + ave1948)/15
A$p8.1972 <- (ave1960*5 + ave1956*4 + ave1952*3 + ave1948*2 + ave1944)/15
A$p8.1968 <- (ave1956*5 + ave1952*4 + ave1948*3 + ave1944*2 + ave1940)/15
A$p8.1964 <- (ave1952*5 + ave1948*4 + ave1944*3 + ave1940*2 + ave1936)/15
A$p8.1960 <- (ave1948*5 + ave1944*4 + ave1940*3 + ave1936*2 + ave1932)/15

A$p12.2008 <- (ave1992*5 + ave1988*4 + ave1984*3 + ave1980*2 + ave1976)/15  ## equivalent to 2024
A$p12.2004 <- (ave1988*5 + ave1984*4 + ave1980*3 + ave1976*2 + ave1972)/15
A$p12.2000 <- (ave1984*5 + ave1980*4 + ave1976*3 + ave1972*2 + ave1968)/15
A$p12.1996 <- (ave1980*5 + ave1976*4 + ave1972*3 + ave1968*2 + ave1964)/15
A$p12.1992 <- (ave1976*5 + ave1972*4 + ave1968*3 + ave1964*2 + ave1960)/15
A$p12.1988 <- (ave1972*5 + ave1968*4 + ave1964*3 + ave1960*2 + ave1956)/15
A$p12.1984 <- (ave1968*5 + ave1964*4 + ave1960*3 + ave1956*2 + ave1952)/15
A$p12.1980 <- (ave1964*5 + ave1960*4 + ave1956*3 + ave1952*2 + ave1948)/15
A$p12.1976 <- (ave1960*5 + ave1956*4 + ave1952*3 + ave1948*2 + ave1944)/15
A$p12.1972 <- (ave1956*5 + ave1952*4 + ave1948*3 + ave1944*2 + ave1940)/15
A$p12.1968 <- (ave1952*5 + ave1948*4 + ave1944*3 + ave1940*2 + ave1936)/15
A$p12.1964 <- (ave1948*5 + ave1944*4 + ave1940*3 + ave1936*2 + ave1932)/15

A$p16.2008 <- (ave1988*5 + ave1984*4 + ave1980*3 + ave1976*2 + ave1972)/15  ## equivalent to 2028
A$p16.2004 <- (ave1984*5 + ave1980*4 + ave1976*3 + ave1972*2 + ave1968)/15
A$p16.2000 <- (ave1980*5 + ave1976*4 + ave1972*3 + ave1968*2 + ave1964)/15
A$p16.1996 <- (ave1976*5 + ave1972*4 + ave1968*3 + ave1964*2 + ave1960)/15
A$p16.1992 <- (ave1972*5 + ave1968*4 + ave1964*3 + ave1960*2 + ave1956)/15
A$p16.1988 <- (ave1968*5 + ave1964*4 + ave1960*3 + ave1956*2 + ave1952)/15
A$p16.1984 <- (ave1964*5 + ave1960*4 + ave1956*3 + ave1952*2 + ave1948)/15
A$p16.1980 <- (ave1960*5 + ave1956*4 + ave1952*3 + ave1948*2 + ave1944)/15
A$p16.1976 <- (ave1956*5 + ave1952*4 + ave1948*3 + ave1944*2 + ave1940)/15
A$p16.1972 <- (ave1952*5 + ave1948*4 + ave1944*3 + ave1940*2 + ave1936)/15
A$p16.1968 <- (ave1948*5 + ave1944*4 + ave1940*3 + ave1936*2 + ave1932)/15

detach(A)




##################################################################
## Find how wrong Burmila's model is on average across the last ##
## 11 elections (1968-2008).    (Old method.)                   ##
##################################################################



attach(A)  ## Matrices of actual and predicted election results
actuals.68 <- cbind(ave2008,ave2004,ave2000,ave1996,ave1992, ave1988, ave1984, ave1980, ave1976, ave1972, ave1968)
predicts.p.68 <- cbind(p.2008,p.2004,p.2000,p.1996,p.1992, p.1988, p.1984, p.1980, p.1976, p.1972, p.1968)
predicts.p4.68 <- cbind(p4.2008,p4.2004,p4.2000,p4.1996,p4.1992, p4.1988, p4.1984, p4.1980, p4.1976, p4.1972, p4.1968)
predicts.p8.68 <- cbind(p8.2008,p8.2004,p8.2000,p8.1996,p8.1992, p8.1988, p8.1984, p8.1980, p8.1976, p8.1972, p8.1968)
predicts.p12.68 <- cbind(p12.2008,p12.2004,p12.2000,p12.1996,p12.1992, p12.1988, p12.1984, p12.1980, p12.1976, p12.1972, p12.1968)
predicts.p16.68 <- cbind(p16.2008,p16.2004,p16.2000,p16.1996,p16.1992, p16.1988, p16.1984, p16.1980, p16.1976, p16.1972, p16.1968)
detach(A)

##################################################################
##              HOW WRONG USING STANDARD DEVIATION              ##
##################################################################



test.wrong <- matrix(NA, nrow=nrow(actuals.68), ncol=5)

for(i in 1:nrow(test.wrong)){
	test.wrong[i,1] <- sqrt(mean((actuals.68[i,] - predicts.p.68[i,])^2, na.rm=T))
	test.wrong[i,2] <- sqrt(mean((actuals.68[i,] - predicts.p4.68[i,])^2, na.rm=T))
	test.wrong[i,3] <- sqrt(mean((actuals.68[i,] - predicts.p8.68[i,])^2, na.rm=T))
	test.wrong[i,4] <- sqrt(mean((actuals.68[i,] - predicts.p12.68[i,])^2, na.rm=T))
	test.wrong[i,5] <- sqrt(mean((actuals.68[i,] - predicts.p16.68[i,])^2, na.rm=T))
	}

mean.diffs.wnv <- c(mean(test.wrong[,1]), mean(test.wrong[,2]), mean(test.wrong[,3]), mean(test.wrong[,4]), mean(test.wrong[,5]))


#           [,1]
#[1,] 0.07460987
#[2,] 0.08188241
#[3,] 0.08327426
#[4,] 0.09046564
#[5,] 0.10289440


##################################################################
##                IOWA and NEW HAMPSHIRE                        ##        
##################################################################

diff.stor <- c()
sds <- seq(1:51)

for(j in 1:nrow(actuals.68)){	
	acts <- actuals.68[j,]		
	preds <- predicts.p.68[j,]	
	
	for(i in 1:1000){
		draw <- sample(1:11, 11, replace=T)
		subseta <- acts[draw]
		subsetb <- preds[draw]
		test <- sqrt(mean((subseta - subsetb)^2, na.rm=T))
		diff.stor[i] <- test
		}
	sds[j] <- mean(diff.stor, na.rm=T)
}

sds

which(A$state=="New Hampshire") #16, 30

attach(A)
NV.p <- (ave2008*5 + ave2004*4 + ave2000*3 + ave1996*2 + ave1992)/15
detach(A)

pdf(file="IA_NH_dists.pdf", width=5, height=5, family="Helvetica", pointsize=10)
plot(density(rnorm(100000, mean=NV.p[16], sd=sds[16])), ylim=c(0,9), col="blue", lty="dotted", lwd=2, main="")
lines(density(rnorm(100000, mean=NV.p[30], sd=sds[30])), col="blue", lty="dashed", lwd=2)
abline(v=.5)
legend("topleft", legend=c("Iowa", "New Hampshire"), col=c("blue", "blue"), lty=c("dotted", "dashed"), lwd=c(2,2))
dev.off()

sds[16]
sds[30]
pnorm(.5, mean=NV.p[16], sd=sds[16])
pnorm(.5, mean=NV.p[30], sd=sds[30])




##################################################################
##        Graph the WNV with 2012 confidence intervals.         ##
##################################################################


unc.graph <- function(sds, NV, mid, name=A$abb){
	dems <- which(NV>=mid)
	reps <- which(NV<mid)
	NV.dems <- NV[dems]
	NV.reps <- NV[reps]
	ub.dem <- NV.dems + (2.36*sds[dems])
	lb.dem <- NV.dems - (2.36*sds[dems])
	ub.rep <- NV.reps + (2.36*sds[reps])
	lb.rep <- NV.reps - (2.36*sds[reps])
	colors.dem <- c()
	colors.rep <- c()
	for(i in 1:length(lb.dem)){
		if(lb.dem[i]<=mid){colors.dem[i]="black"}
		else{colors.dem[i]="blue"}
		}
	for(i in 1:length(lb.rep)){
		if(ub.rep[i]>mid){colors.rep[i]="black"}
		else{colors.rep[i]="red"}
		}
	plot(x=NV.dems, y=dems, pch=20, col="blue", ylab="States, Ordered Alphabetically", xlim=c(mid-.5,mid+.5), ylim=c(0,52), xlab="Democratic Vote Share", yaxt="n")
	text(lb.dem-.02, dems, labels=name[dems], cex=.7, col=colors.dem)
	points(x=NV.reps, y=reps, pch=20, col="red")
	text(lb.rep-.02, reps, labels=name[reps], cex=.7, col=colors.rep)
	segments(x0=lb.dem, x1=ub.dem, y0=dems, y1=dems, col="blue")
	segments(x0=lb.rep, x1=ub.rep, y0=reps, y1=reps, col="red")
	abline(v=mid, lty="dotted")
	}


unc.graph(sds, NV.p, .5)

##################################################################
## Calculate relative state-level vote.                         ##
##################################################################

attach(A)  ## US-wide vote share.
us2008 <- sum(d2008)/(sum(d2008+r2008))
us2004 <- sum(d2004)/(sum(d2004+r2004))
us2000 <- sum(d2000)/(sum(d2000+r2000))
us1996 <- sum(d1996)/(sum(d1996+r1996))
us1992 <- sum(d1992)/(sum(d1992+r1992))
us1988 <- sum(d1988)/(sum(d1988+r1988))
us1984 <- sum(d1984)/(sum(d1984+r1984))
us1980 <- sum(d1980)/(sum(d1980+r1980))
us1976 <- sum(d1976)/(sum(d1976+r1976))
us1972 <- sum(d1972)/(sum(d1972+r1972))
us1968 <- sum(d1968)/(sum(d1968+r1968))
us1964 <- sum(d1964)/(sum(d1964+r1964))
us1960 <- sum(d1960, na.rm=T)/(sum(d1960, na.rm=T)+sum(r1960, na.rm=T))
us1956 <- sum(d1956, na.rm=T)/(sum(d1956, na.rm=T)+sum(r1956, na.rm=T))
us1952 <- sum(d1952, na.rm=T)/(sum(d1952, na.rm=T)+sum(r1952, na.rm=T))
us1948 <- sum(d1948, na.rm=T)/(sum(d1948, na.rm=T)+sum(r1948, na.rm=T))
us1944 <- sum(d1944, na.rm=T)/(sum(d1944, na.rm=T)+sum(r1944, na.rm=T))
us1940 <- sum(d1940, na.rm=T)/(sum(d1940, na.rm=T)+sum(r1940, na.rm=T))
us1936 <- sum(d1936, na.rm=T)/(sum(d1936, na.rm=T)+sum(r1936, na.rm=T))
us1932 <- sum(d1932, na.rm=T)/(sum(d1932, na.rm=T)+sum(r1932, na.rm=T))
detach(A)

usave <- c(us2008, us2004, us2000, us1996, us1992, us1988, us1984, us1980, us1976, us1972, us1968, us1964, us1960, us1956, us1952, us1948, us1944, us1940, us1936, us1932)

mean(usave); sd(usave)

diffs <- data.frame()
for(i in 1:nrow(A)){
	for(j in 1:length(usave)){
		diffs[i,j] <- A[i,j+42]-usave[j]
		}	
	}

colnames(diffs) <- c("diff2008", "diff2004", "diff2000", "diff1996", "diff1992", "diff1988", "diff1984", "diff1980", "diff1976", "diff1972", "diff1968", "diff1964", "diff1960", "diff1956", "diff1952", "diff1948", "diff1944", "diff1940", "diff1936", "diff1932")


##################################################################
##            Make predictions based on relative vote           ##
##################################################################

attach(diffs)
diffs$prw.2008 <- (diff2004*5 + diff2000*4 + diff1996*3 + diff1992*2 + diff1988)/15
diffs$prw.2004 <- (diff2000*5 + diff1996*4 + diff1992*3 + diff1988*2 + diff1984)/15
diffs$prw.2000 <- (diff1996*5 + diff1992*4 + diff1988*3 + diff1984*2 + diff1980)/15
diffs$prw.1996 <- (diff1992*5 + diff1988*4 + diff1984*3 + diff1980*2 + diff1976)/15
diffs$prw.1992 <- (diff1988*5 + diff1984*4 + diff1980*3 + diff1976*2 + diff1972)/15
diffs$prw.1988 <- (diff1984*5 + diff1980*4 + diff1976*3 + diff1972*2 + diff1968)/15
diffs$prw.1984 <- (diff1980*5 + diff1976*4 + diff1972*3 + diff1968*2 + diff1964)/15
diffs$prw.1980 <- (diff1976*5 + diff1972*4 + diff1968*3 + diff1964*2 + diff1960)/15
diffs$prw.1976 <- (diff1972*5 + diff1968*4 + diff1964*3 + diff1960*2 + diff1956)/15
diffs$prw.1972 <- (diff1968*5 + diff1964*4 + diff1960*3 + diff1956*2 + diff1952)/15
diffs$prw.1968 <- (diff1964*5 + diff1960*4 + diff1956*3 + diff1952*2 + diff1948)/15
diffs$prw.1964 <- (diff1960*5 + diff1956*4 + diff1952*3 + diff1948*2 + diff1944)/15
diffs$prw.1960 <- (diff1956*5 + diff1952*4 + diff1948*3 + diff1944*2 + diff1940)/15
diffs$prw.1956 <- (diff1952*5 + diff1948*4 + diff1944*3 + diff1940*2 + diff1936)/15
diffs$prw.1952 <- (diff1948*5 + diff1944*4 + diff1940*3 + diff1936*2 + diff1932)/15

diffs$prw4.2008 <- (diff2000*5 + diff1996*4 + diff1992*3 + diff1988*2 + diff1984)/15
diffs$prw4.2004 <- (diff1996*5 + diff1992*4 + diff1988*3 + diff1984*2 + diff1980)/15
diffs$prw4.2000 <- (diff1992*5 + diff1988*4 + diff1984*3 + diff1980*2 + diff1976)/15
diffs$prw4.1996 <- (diff1988*5 + diff1984*4 + diff1980*3 + diff1976*2 + diff1972)/15
diffs$prw4.1992 <- (diff1984*5 + diff1980*4 + diff1976*3 + diff1972*2 + diff1968)/15
diffs$prw4.1988 <- (diff1980*5 + diff1976*4 + diff1972*3 + diff1968*2 + diff1964)/15
diffs$prw4.1984 <- (diff1976*5 + diff1972*4 + diff1968*3 + diff1964*2 + diff1960)/15
diffs$prw4.1980 <- (diff1972*5 + diff1968*4 + diff1964*3 + diff1960*2 + diff1956)/15
diffs$prw4.1976 <- (diff1968*5 + diff1964*4 + diff1960*3 + diff1956*2 + diff1952)/15
diffs$prw4.1972 <- (diff1964*5 + diff1960*4 + diff1956*3 + diff1952*2 + diff1948)/15
diffs$prw4.1968 <- (diff1960*5 + diff1956*4 + diff1952*3 + diff1948*2 + diff1944)/15
diffs$prw4.1964 <- (diff1956*5 + diff1952*4 + diff1948*3 + diff1944*2 + diff1940)/15
diffs$prw4.1960 <- (diff1952*5 + diff1948*4 + diff1944*3 + diff1940*2 + diff1936)/15
diffs$prw4.1956 <- (diff1948*5 + diff1944*4 + diff1940*3 + diff1936*2 + diff1932)/15

diffs$prw8.2008 <- (diff1996*5 + diff1992*4 + diff1988*3 + diff1984*2 + diff1980)/15
diffs$prw8.2004 <- (diff1992*5 + diff1988*4 + diff1984*3 + diff1980*2 + diff1976)/15
diffs$prw8.2000 <- (diff1988*5 + diff1984*4 + diff1980*3 + diff1976*2 + diff1972)/15
diffs$prw8.1996 <- (diff1984*5 + diff1980*4 + diff1976*3 + diff1972*2 + diff1968)/15
diffs$prw8.1992 <- (diff1980*5 + diff1976*4 + diff1972*3 + diff1968*2 + diff1964)/15
diffs$prw8.1988 <- (diff1976*5 + diff1972*4 + diff1968*3 + diff1964*2 + diff1960)/15
diffs$prw8.1984 <- (diff1972*5 + diff1968*4 + diff1964*3 + diff1960*2 + diff1956)/15
diffs$prw8.1980 <- (diff1968*5 + diff1964*4 + diff1960*3 + diff1956*2 + diff1952)/15
diffs$prw8.1976 <- (diff1964*5 + diff1960*4 + diff1956*3 + diff1952*2 + diff1948)/15
diffs$prw8.1972 <- (diff1960*5 + diff1956*4 + diff1952*3 + diff1948*2 + diff1944)/15
diffs$prw8.1968 <- (diff1956*5 + diff1952*4 + diff1948*3 + diff1944*2 + diff1940)/15
diffs$prw8.1964 <- (diff1952*5 + diff1948*4 + diff1944*3 + diff1940*2 + diff1936)/15
diffs$prw8.1960 <- (diff1948*5 + diff1944*4 + diff1940*3 + diff1936*2 + diff1932)/15

diffs$prw12.2008 <- (diff1992*5 + diff1988*4 + diff1984*3 + diff1980*2 + diff1976)/15
diffs$prw12.2004 <- (diff1988*5 + diff1984*4 + diff1980*3 + diff1976*2 + diff1972)/15
diffs$prw12.2000 <- (diff1984*5 + diff1980*4 + diff1976*3 + diff1972*2 + diff1968)/15
diffs$prw12.1996 <- (diff1980*5 + diff1976*4 + diff1972*3 + diff1968*2 + diff1964)/15
diffs$prw12.1992 <- (diff1976*5 + diff1972*4 + diff1968*3 + diff1964*2 + diff1960)/15
diffs$prw12.1988 <- (diff1972*5 + diff1968*4 + diff1964*3 + diff1960*2 + diff1956)/15
diffs$prw12.1984 <- (diff1968*5 + diff1964*4 + diff1960*3 + diff1956*2 + diff1952)/15
diffs$prw12.1980 <- (diff1964*5 + diff1960*4 + diff1956*3 + diff1952*2 + diff1948)/15
diffs$prw12.1976 <- (diff1960*5 + diff1956*4 + diff1952*3 + diff1948*2 + diff1944)/15
diffs$prw12.1972 <- (diff1956*5 + diff1952*4 + diff1948*3 + diff1944*2 + diff1940)/15
diffs$prw12.1968 <- (diff1952*5 + diff1948*4 + diff1944*3 + diff1940*2 + diff1936)/15
diffs$prw12.1964 <- (diff1948*5 + diff1944*4 + diff1940*3 + diff1936*2+ diff1932)/15

diffs$prw16.2008 <- (diff1988*5 + diff1984*4 + diff1980*3 + diff1976*2 + diff1972)/15
diffs$prw16.2004 <- (diff1984*5 + diff1980*4 + diff1976*3 + diff1972*2 + diff1968)/15
diffs$prw16.2000 <- (diff1980*5 + diff1976*4 + diff1972*3 + diff1968*2 + diff1964)/15
diffs$prw16.1996 <- (diff1976*5 + diff1972*4 + diff1968*3 + diff1964*2 + diff1960)/15
diffs$prw16.1992 <- (diff1972*5 + diff1968*4 + diff1964*3 + diff1960*2 + diff1956)/15
diffs$prw16.1988 <- (diff1968*5 + diff1964*4 + diff1960*3 + diff1956*2 + diff1952)/15
diffs$prw16.1984 <- (diff1964*5 + diff1960*4 + diff1956*3 + diff1952*2 + diff1948)/15
diffs$prw16.1980 <- (diff1960*5 + diff1956*4 + diff1952*3 + diff1948*2 + diff1944)/15
diffs$prw16.1976 <- (diff1956*5 + diff1952*4 + diff1948*3 + diff1944*2 + diff1940)/15
diffs$prw16.1972 <- (diff1952*5 + diff1948*4 + diff1944*3 + diff1940*2 + diff1936)/15
diffs$prw16.1968 <- (diff1948*5 + diff1944*4 + diff1940*3 + diff1936*2 + diff1932)/15
detach(diffs)

## Matrices of actual and predicted Dem vote deviation

attach(diffs)  
actuals.r68 <- cbind(diff2008,diff2004,diff2000,diff1996,diff1992, diff1988, diff1984, diff1980, diff1976, diff1972, diff1968)
predicts.prw.68 <- cbind(prw.2008,prw.2004,prw.2000,prw.1996,prw.1992, prw.1988, prw.1984, prw.1980, prw.1976, prw.1972, prw.1968)
predicts.prw4.68 <- cbind(prw4.2008,prw4.2004,prw4.2000,prw4.1996,prw4.1992, prw4.1988, prw4.1984, prw4.1980, prw4.1976, prw4.1972, prw4.1968)
predicts.prw8.68 <- cbind(prw8.2008,prw8.2004,prw8.2000,prw8.1996,prw8.1992, prw8.1988, prw8.1984, prw8.1980, prw8.1976, prw8.1972, prw8.1968)
predicts.prw12.68 <- cbind(prw12.2008,prw12.2004,prw12.2000,prw12.1996,prw12.1992, prw12.1988, prw12.1984, prw12.1980, prw12.1976, prw12.1972, prw12.1968)
predicts.prw16.68 <- cbind(prw16.2008,prw16.2004,prw16.2000,prw16.1996,prw16.1992, prw16.1988, prw16.1984, prw16.1980, prw16.1976, prw16.1972, prw16.1968)
detach(diffs)


test.wrong <- matrix(NA, nrow=nrow(actuals.68), ncol=5)

for(i in 1:nrow(test.wrong)){
	test.wrong[i,1] <- sqrt(mean((actuals.r68[i,] - predicts.prw.68[i,])^2, na.rm=T))
	test.wrong[i,2] <- sqrt(mean((actuals.r68[i,] - predicts.prw4.68[i,])^2, na.rm=T))
	test.wrong[i,3] <- sqrt(mean((actuals.r68[i,] - predicts.prw8.68[i,])^2, na.rm=T))
	test.wrong[i,4] <- sqrt(mean((actuals.r68[i,] - predicts.prw12.68[i,])^2, na.rm=T))
	test.wrong[i,5] <- sqrt(mean((actuals.r68[i,] - predicts.prw16.68[i,])^2, na.rm=T))
	}

mean.diffs.dev <- c(mean(test.wrong[,1]), mean(test.wrong[,2]), mean(test.wrong[,3]), mean(test.wrong[,4]), mean(test.wrong[,5]))
storage <- cbind(mean.diffs.wnv, mean.diffs.dev)
	
rownames(storage) <- c("One Election Ahead", "Two Elections Ahead", "Three Elections Ahead", "Four Elections Ahead", "Five Elections Ahead")

xtable(storage, digits=4)

##################################################################
## Calculating electoral votes - using Census predictions and   ##
## ACS estimates.                                               ##
##################################################################

multiplier <- matrix(NA, 59, 1)	 
for(i in 1:59){
	out<-1/sqrt((i+1)*i) 
	multiplier[i,] <- out 
	}

ev.function <- function(x){		
	melt <- melt(x)
	order <- melt[order(melt$value, decreasing=TRUE),]
	final <- order[1:385,1]
	out <- table(final)+3		 
	return(out)}				

m.2000 <- data.frame(A$state)
for(i in 1:59){
	m.2000[,paste(sep="","var",i)] <- B$evp2000*multiplier[i,]
	}
m.2010 <- data.frame(A$state)
for(i in 1:59){
	m.2010[,paste(sep="","var",i)] <- B$p2010*multiplier[i,]
	}
m.2009 <- data.frame(A$state)
for(i in 1:59){
	m.2009[,paste(sep="","var",i)] <- B$e2009*multiplier[i,]
	}

B$ev2000 <- ev.function(m.2000)	
B$ev2010 <- ev.function(m.2010)	
B$ev2009 <- ev.function(m.2009)

B$ev2010-B$ev2000
B$ev2009-B$ev2000

##################################################################
## Using the PRW and its 2012 standard error, simulate the net  ##
## gain or loss in electoral votes for the Democratic party.    ##        
##################################################################


attach(diffs)
prw <- (diff2008*5 + diff2004*4 + diff2000*3 + diff1996*2 + diff1992)/15
detach(diffs)


diff.stor <- seq(1:51)
sds <- c()


for(j in 1:nrow(actuals.r68)){	
	acts <- actuals.r68[j,]		
	preds <- predicts.prw.68[j,]	
	
	for(i in 1:1000){
		draw <- sample(1:11, 11, replace=T)
		subseta <- acts[draw]
		subsetb <- preds[draw]
		test <- sqrt(mean((subseta - subsetb)^2, na.rm=T))
		diff.stor[i] <- test
		}
	sds[j] <- mean(diff.stor, na.rm=T)
}

mean(sds)

Dgainloss.p <- c()

for(i in 1:10000){
	out <- rep(0, 15)
	draw.USA <- rnorm(1, mean(usave), sd(usave))
	draw.AL <- rnorm(1, mean=prw[1], sds[1]) + draw.USA
		if(draw.AL>=.5){out[1]<-(-1)}
		else{out[1]<-(1)}
	draw.AZ <- rnorm(1, mean=prw[3], sds[3]) + draw.USA
		if(draw.AZ>=.5){out[2]<-1}
		else{out[2]<-(-1)}
	draw.CA <- rnorm(1, mean=prw[5], sds[5]) + draw.USA
		if(draw.CA>=.5){out[3]<-1}
		else{out[3]<-(-1)}
	draw.FL <- rnorm(1, mean=prw[10], sds[10]) + draw.USA
		if(draw.FL>=.5){out[4]<-2}
		else{out[4]<-(-2)}
	draw.GA <- rnorm(1, mean=prw[11], sds[11]) + draw.USA
		if(draw.GA>=.5){out[5]<-1}
		else{out[5]<-(-1)}
	draw.IL <- rnorm(1, mean=prw[14], sds[14]) + draw.USA
		if(draw.IL>=.5){out[6]<-(-1)}
		else{out[6]<-(1)}
	draw.IA <- rnorm(1, mean=prw[16], sds[16]) + draw.USA
		if(draw.IA>=.5){out[7]<-(-1)}
		else{out[7]<-(1)}
	draw.MA <- rnorm(1, mean=prw[22], sds[22]) + draw.USA
		if(draw.MA>=.5){out[8]<-(-1)}
		else{out[8]<-(1)}
	draw.MO <- rnorm(1, mean=prw[26], sds[26]) + draw.USA
		if(draw.MO>=.5){out[9]<-(-1)}
		else{out[9]<-(1)}
	draw.NV <- rnorm(1, mean=prw[29], sds[29]) + draw.USA
		if(draw.NV>=.5){out[10]<-1}
		else{out[10]<-(-1)}
	draw.NY <- rnorm(1, mean=prw[33], sds[33]) + draw.USA
		if(draw.NY>=.5){out[11]<-(-2)}
		else{out[11]<-(2)}
	draw.OH <- rnorm(1, mean=prw[36], sds[36]) + draw.USA
		if(draw.OH>=.5){out[12]<-(-2)}
		else{out[12]<-(2)}
	draw.PA <- rnorm(1, mean=prw[39], sds[39]) + draw.USA
		if(draw.PA>=.5){out[13]<-(-1)}
		else{out[13]<-(1)}
	draw.TX <- rnorm(1, mean=prw[44], sds[44]) + draw.USA
		if(draw.TX>=.5){out[14]<-3}
		else{out[14]<-(-3)}
	draw.UT <- rnorm(1, mean=prw[45], sds[45]) + draw.USA
		if(draw.UT>=.5){out[15]<-1}
		else{out[15]<-(-1)}
	Dgainloss.p[i] <- sum(out)
	}

a <- table(Dgainloss.p)
b <- a; b[10] <- 0; b[11] <- 0; b[12] <- 0
c <- rep(0,12); c[10] <- a[10]
d <- c(rep(0,10), a[11:12])
a <- rbind(b,c, d)

pdf(file="Dgainloss.pdf", width=5, height=5, family="Helvetica", pointsize=10)
barplot(a, beside=FALSE, col=c("brown1", "white", "cornflowerblue"))
dev.off()

mean(Dgainloss.p)
length(which(Dgainloss.p>=0))/5000
length(which(Dgainloss.p<(-9)))/5000


##################################################################
## 1988 Maps - Burmila and PRW Categories                       ##
##################################################################

cat.colors <- function(x){		
	if(x<.455){					
		out <- "red"
		}
	else if(x>.455 & x<.495){  
		out <- "brown1" 
		}  
	else if(x>.495 & x<.515){
		out <- "white"
		}
	else if(x>.515 & x<.5501){  
		out <- "cornflowerblue"
		}
	else{
		out <- "blue"
		}
	return(out)
	}

cat.colors.r <- function(x){
	if(x<(-.10)){out <- "red"}
	else if(x>-.10 & x<(-.02)){out <- "brown1"}
	else if(x>-.02 & x<.02){out <- "white"}
	else if(x>.02 & x<.10){out <- "cornflowerblue"}
	else{out <- "blue"}
	return(out)
	}

map.1992 <- matrix(NA, nrow=51, ncol=3)
for(i in 1:nrow(map.1992)){
	map.1992[i,1] <- cat.colors(A$p.1992[i])
	map.1992[i,2] <- cat.colors.r(diffs$prw.1992[i])
	}

map.1992[2:10,] <- map.1992[3:11,]
map.1992[11:21,] <- map.1992[13:23,]
map.1992[21,] <- map.1992[20,]
map.1992 <- rbind(map.1992[1:23,], map.1992[23,], map.1992[24:51,]) 
map.1992 <- rbind(map.1992[1:34,], map.1992[34,], map.1992[34,], map.1992[34,], map.1992[35:52,])
map.1992 <- rbind(map.1992[1:38,], map.1992[38,], map.1992[38,], map.1992[39:55,])
map.1992 <- rbind(map.1992[1:53,], map.1992[53,], map.1992[53,], map.1992[54:57,])
map.1992 <- rbind(map.1992[1:56,], map.1992[56,], map.1992[56,], map.1992[56,], map.1992[56,], map.1992[57:59,])

pdf(file="WNV1988.pdf", width=5, height=5, family="Helvetica", pointsize=10)
map("state", fill=TRUE, col=map.1992[,1])
legend("bottomright", fill=c("blue", "cornflowerblue", "white", "brown1", "red"), legend=c("SD", "LD", "TU", "LR", "SR"), cex=.6)
dev.off()

pdf(file="WD_1988.pdf", width=5, height=5, family="Helvetica", pointsize=10)
map("state", fill=TRUE, col=map.1992[,2])
legend("bottomright", fill=c("blue", "cornflowerblue", "white", "brown1", "red"), legend=c("SD", "LD", "TU", "LR", "SR"), cex=.6)
dev.off()

category <- function(x){		
	if(x<.455){out <- "SR"}
	else if(x>.455 & x<.495){out <- "LR"}  
	else if(x>.495 & x<.515){out <- "TU"}
	else if(x>.515 & x<.5501){out <- "LD"}
	else{out <- "SD"}
	return(out)
	}

category.r <- function(x){
	if(x<(-.10)){out <- "SR"}
	else if(x>-.10 & x<(-.02)){out <- "LR"}
	else if(x>-.02 & x<.02){out <- "TU"}
	else if(x>.02 & x<.10){out <- "LD"}
	else{out <- "SD"}
	return(out)
	}

cat.1992 <- data.frame()
for(i in 1:nrow(A)){
	cat.1992[i,1] <- category(A$p.1992[i])
	cat.1992[i,2] <- category.r(diffs$prw.1992[i])
	cat.1992[i,3] <- category(A$ave1992[i])
	cat.1992[i,4] <- category(A$ave1996[i])
	cat.1992[i,5] <- category(A$ave2000[i])
	cat.1992[i,6] <- category(A$ave2004[i])
	cat.1992[i,7] <- category(A$ave2008[i])
	}
cat.1992 <- cbind(cat.1992, B$ev1990, B$ev2000)

pred <- with(cat.1992, tapply(B$ev1990,V1,sum))
pred.r <- with(cat.1992, tapply(B$ev1990, V2, sum))
a1992 <- with(cat.1992, tapply(B$ev1990,V3,sum))
a1996 <- with(cat.1992, tapply(B$ev1990,V4,sum))
a2000 <- with(cat.1992, tapply(B$ev1990,V5,sum))
a2004 <- with(cat.1992, tapply(Freq,V6,sum))
a2008 <- with(cat.1992, tapply(Freq,V7,sum))

a1992-pred
a1996-pred
a2000-pred
a2004-pred
a2008-pred

a1992-pred.r
a1996-pred.r
a2000-pred.r
a2004-pred.r
a2008-pred.r

##################################################################
## Comparing Census 1990 predictions to actual EV distributions ##
##################################################################

m.2000 <- data.frame(A$state)
for(i in 1:59){
	m.2000[,paste(sep="","var",i)] <- B$pe2000*multiplier[i,]
	}
m.2010 <- data.frame(A$state)
for(i in 1:59){
	m.2010[,paste(sep="","var",i)] <- B$pe2010*multiplier[i,]
	}

B$ev2000.p <- ev.function(m.2000)
B$error.2000 <- B$ev2000-B$ev2000.p ## 6 states wrong
B$ev2010.p <- ev.function(m.2010)
B$error.2010 <- B$ev2009-B$ev2010.p ## 14 states wrong

##################################################################
## Electoral vote test                                          ##
##################################################################

attach(B)
ev.mat <- cbind(ev1930, ev1940, ev1950, ev1960, ev1970, ev1980, ev1990, ev2000)
detach(B)

attach(A)
vote.mat <- cbind(ave1932, ave1936, ave1940, ave1944, ave1948, ave1952, ave1956, ave1960, ave1964, ave1968, ave1972, ave1976, ave1980, ave1984, ave1988, ave1992, ave1996, ave2000, ave2004, ave2008)
detach(A)

dem <- matrix(0, nrow(vote.mat), nrow(vote.mat))
for(i in 1:nrow(vote.mat)){
	for(j in 1:ncol(vote.mat)){
		if(is.na(vote.mat[i,j])){dem[i,j]=0}
		else if(vote.mat[i,j]>.5){dem[i,j]=1}
		}
	}
ev.sim <- matrix(NA, ncol(ev.mat), ncol(vote.mat))
for(i in 1:ncol(vote.mat)){
	a <- dem[,i]
	b <- c()
	for(j in 1:ncol(ev.mat)){
		b[j] <- sum(a* ev.mat[,j], na.rm=T)
		}
	ev.sim[,i] <- b
	}

colnames(ev.sim) <- c("D1932", "D1936", "D1940", "D1944", "D1948", "D1952", "D1956", "D1960", "D1964", "D1968", "D1972", "D1976", "D1980", "D1984", "D1988", "D1992", "D1996", "D2000", "D2004", "D2008")

## Elections where EV matters: 1960, 2000, 2004 (15% of elections)
## In 1960, R would have won with any distribution after 1960 (5 of 8 distributions)
## In 2000, D would have won with 4 of 8 distributions (all elections before 2000, if you add in AK, DC, HI)
## In 2004, D would have won with 2 of 8 distributions
## Only other particularly close one is 1976

##################################################################
## Demographic Model using TX County Information                ##
##################################################################

pop2008 <- read.csv("alldata.csv")
# Estimates of population by age, race, gender, by county.  2008 only.  http://txsdc.utsa.edu/tpepp/txpopest.php
vote2008 <- read.csv("2008vote_bycounty.csv")
vote2004 <- read.csv("2004vote_bycounty.csv")
vote2004 <- vote2004[1:254,]
vote2000 <- read.csv("2000vote_bycounty.csv")
vote1996 <- read.csv("1996vote_bycounty.csv")
vote1992 <- read.csv("1992vote_bycounty.csv")
## http://elections.sos.state.tx.us/elchist.exe
acs <- read.csv("Data_InclMedInc.csv")

###################   2008 ESTIMATES   ###################

aa.pop2008 <- pop2008[pop2008$Age=="All Ages",]

perc.hisp <- aa.pop2008[1,16]/aa.pop2008[1,4]
toth.2008 <- aa.pop2008[1,16]
tot.2008 <- aa.pop2008[1,4]
perchisp2008 <- toth.2008/tot.2008

colnames(aa.pop2008)

### adding 2008 estimates to the dataset 

aa.pop2008 <- aa.pop2008[-1,]
aa.pop2008 <- aa.pop2008[,c(4,7,10,13,16)]
colnames(aa.pop2008) <- c("aa2008", "aaw2008", "aab2008", "aao2008", "aah2008")
head(aa.pop2008)
dim(aa.pop2008)

hisp2008 <- aa.pop2008$aah2008
tot2008 <- aa.pop2008$aa2008
ph2008 <- hisp2008/tot2008

################### COUNTY MEDIAN INCOME ###################

acs.tx <- acs[acs$ANSI>=48000,]
acs.tx <- acs.tx[acs.tx$ANSI<49000,]
dim(acs.tx)
medianinc <- acs.tx$IPE010207[-1]


################### CALCULATE COUNTRY DIFFS ###################

dvotesharetx <- sum(vote2008$Dvote)/sum(vote2008$Dvote + vote2008$Rvote) 

dvote2008 <- vote2008[,7]
diff.txcty08 <- dvote2008-dvotesharetx

dvoteshare04 <- sum(vote2004$DEM)/sum(vote2004$DEM+vote2004$REP)
dvotesharecty <- vote2004$DEM/(vote2004$DEM+vote2004$REP)
diff.txcty04 <- dvotesharecty-dvoteshare04

dvoteshare00 <- sum(vote2000$DEM)/sum(vote2000$DEM+vote2000$REP)
dvotesharecty <- vote2000$DEM/(vote2000$DEM+vote2000$REP)
diff.txcty00 <- dvotesharecty-dvoteshare00

dvoteshare96 <- sum(vote1996$DEM)/sum(vote1996$DEM+vote1996$REP)
dvotesharecty <- vote1996$DEM/(vote1996$DEM+vote1996$REP)
diff.txcty96 <- dvotesharecty-dvoteshare96

dvoteshare92 <- sum(vote1992$DEM)/sum(vote1992$DEM+vote1992$REP)
dvotesharecty <- vote1992$DEM/(vote1992$DEM+vote1992$REP)
diff.txcty92 <- dvotesharecty-dvoteshare92


################### CALCULATE DEMOG DATA ###################

perc.hisp.2008 <- aa.pop2008$aah2008/aa.pop2008$aa2008

perc.b.2008 <- aa.pop2008$aab2008/aa.pop2008$aa2008

pop2008.young <- subset(pop2008, pop2008$Age=="18 Years" | pop2008$Age=="19 Years" | pop2008$Age=="20 Years" | pop2008$Age=="21 Years" | pop2008$Age=="22 Years"| pop2008$Age=="23 Years"| pop2008$Age=="24 Years"| pop2008$Age=="25 Years"| pop2008$Age=="26 Years"| pop2008$Age=="27 Years"| pop2008$Age=="28 Years"| pop2008$Age=="29 Years"| pop2008$Age=="30 Years")

pop2008.seniors <- subset(pop2008, pop2008$Age=="65 Years" | pop2008$Age=="66 Years" | pop2008$Age=="67 Years" | pop2008$Age=="68 Years" | pop2008$Age=="69 Years" | pop2008$Age=="70 Years" | pop2008$Age=="71 Years" | pop2008$Age=="72 Years" | pop2008$Age=="73 Years" | pop2008$Age=="74 Years" | pop2008$Age=="75 Years" | pop2008$Age=="76 Years" | pop2008$Age=="77 Years" | pop2008$Age=="78 Years" | pop2008$Age=="79 Years" | pop2008$Age=="80 Years" | pop2008$Age=="81 Years" | pop2008$Age=="82 Years" | pop2008$Age=="83 Years" | pop2008$Age=="84 Years" | pop2008$Age=="85 Years +")


COUNTYNAME <- unique(pop2008.seniors$County)
colnames(pop2008.seniors)

sumytoo <- c()

test <- function(x,y){
	for(i in 1:length(x)){
		ytoo <- subset(y, y[,1]==x[i])
		sumytoo[i] <- sum(ytoo[,4])
		}
		return(sumytoo)
	}

seniors <- test(COUNTYNAME, pop2008.seniors)
young <- test(COUNTYNAME, pop2008.young)
perc.seniors <- seniors[-1]/aa.pop2008$aa2008
perc.young <- young[-1]/aa.pop2008$aa2008


################### REGRESSIONS ###################

tx.2008.mat <- cbind(diff.txcty08, diff.txcty04, diff.txcty00, diff.txcty96, diff.txcty92, medianinc, perc.hisp.2008, perc.b.2008, perc.young, perc.seniors)

xtable(summary(lm(diff.txcty08 ~ diff.txcty92 + perc.hisp.2008 + perc.b.2008+ perc.seniors + perc.young)))

xtable(summary(lm(diff.txcty08 ~ diff.txcty92 + medianinc + perc.hisp.2008 + perc.b.2008+ perc.seniors + perc.young)))

xtable(summary(lm(diff.txcty08 ~ diff.txcty04 + perc.hisp.2008 + perc.b.2008+ perc.seniors + perc.young)))

xtable(summary(lm(diff.txcty08 ~ diff.txcty04 + medianinc + perc.hisp.2008 + perc.b.2008+ perc.seniors + perc.young)))

summary(lm(diff.txcty08 ~ diff.txcty96 + medianinc + perc.hisp.2008 + perc.b.2008+ perc.seniors + perc.young))

summary(lm(diff.txcty08 ~ diff.txcty00 + medianinc + perc.hisp.2008 + perc.b.2008+ perc.seniors + perc.young))

##################################################################
## Cut-off for reach of model                                   ##       
##################################################################

sds4 <- c()
for(j in 1:nrow(actuals.68)){	
	acts <- actuals.r68[j,]		
	preds <- predicts.prw4.68[j,]	
	
	for(i in 1:1000){
		draw <- sample(1:11, 11, replace=T)
		subseta <- acts[draw]
		subsetb <- preds[draw]
		test <- sqrt(mean((subseta - subsetb)^2, na.rm=T))
		diff.stor[i] <- test
		}
	sds4[j] <- mean(diff.stor, na.rm=T)
}

sds8 <- c()
for(j in 1:nrow(actuals.68)){	
	acts <- actuals.r68[j,]		
	preds <- predicts.prw8.68[j,]	
	
	for(i in 1:1000){
		draw <- sample(1:11, 11, replace=T)
		subseta <- acts[draw]
		subsetb <- preds[draw]
		test <- sqrt(mean((subseta - subsetb)^2, na.rm=T))
		diff.stor[i] <- test
		}
	sds8[j] <- mean(diff.stor, na.rm=T)
}

sds12 <- c()
for(j in 1:nrow(actuals.68)){	
	acts <- actuals.r68[j,]		
	preds <- predicts.prw12.68[j,]	
	
	for(i in 1:1000){
		draw <- sample(1:11, 11, replace=T)
		subseta <- acts[draw]
		subsetb <- preds[draw]
		test <- sqrt(mean((subseta - subsetb)^2, na.rm=T))
		diff.stor[i] <- test
		}
	sds12[j] <- mean(diff.stor, na.rm=T)
}

sds16 <- c()
for(j in 1:nrow(actuals.68)){	
	acts <- actuals.r68[j,]		
	preds <- predicts.prw16.68[j,]	
	
	for(i in 1:1000){
		draw <- sample(1:11, 11, replace=T)
		subseta <- acts[draw]
		subsetb <- preds[draw]
		test <- sqrt(mean((subseta - subsetb)^2, na.rm=T))
		diff.stor[i] <- test
		}
	sds16[j] <- mean(diff.stor, na.rm=T)
}

pct <- c(); pct.4 <- c(); pct.8 <- c(); pct.12 <- c(); pct.16 <- c()

for(i in 1:nrow(A)){
	pct[i] <- pnorm(0, mean=prw[i], sd=sds[i])
	pct.4[i] <- pnorm(0, mean=prw[i], sd=sds4[i])
	pct.8[i] <- pnorm(0, mean=prw[i], sd=sds8[i])
	pct.12[i] <- pnorm(0, mean=prw[i], sd=sds12[i])
	pct.16[i] <- pnorm(0, mean=prw[i], sd=sds16[i])	}
	
pct <- data.frame(A$state, round(pct,3), round(pct.4,3), round(pct.8,3), round(pct.12,3), round(pct.16,3))
colnames(pct) <- c("state", "pct", "pct.4", "pct.8", "pct.12", "pct.16")

mean(pct[,2]); mean(pct[,3]); mean(pct[,4]); mean(pct[,5]); mean(pct[,6])

length(pct[which(pct[,2]>.95 | pct[,2]<.05),2])
length(pct[which(pct[,3]>.95 | pct[,3]<.05),3])
length(pct[which(pct[,4]>.95 | pct[,4]<.05),4])
length(pct[which(pct[,5]>.95 | pct[,5]<.05),5])
length(pct[which(pct[,6]>.95 | pct[,6]<.05),6])


